!> \file mynnedmf_wrapper.F90
!!  This file contains all of the code related to running the MYNN 
!! eddy-diffusivity mass-flux scheme. 

!> The following references best describe the code within
!!    Olson et al. (2019, NOAA Technical Memorandum)
!!    Nakanishi and Niino (2009) \cite NAKANISHI_2009
      MODULE mynnedmf_wrapper

      contains

!> \section arg_table_mynnedmf_wrapper_init Argument Table
!! \htmlinclude mynnedmf_wrapper_init.html
!!
      subroutine mynnedmf_wrapper_init (          &
        &  con_cp, con_grav, con_rd, con_rv,      &
        &  con_cpv, con_cliq, con_cice, con_rcp,  &
        &  con_XLV, con_XLF, con_p608, con_ep2,   &
        &  con_karman, con_t0c,                   &
        &  do_mynnedmf, lheatstrg,                &
        &  errmsg, errflg                         )

        use machine,  only : kind_phys
        use bl_mynn_common

        implicit none
        
        logical,             intent(in)  :: do_mynnedmf
        logical,             intent(in)  :: lheatstrg
        character(len=*),    intent(out) :: errmsg
        integer,             intent(out) :: errflg

        real(kind=kind_phys),intent(in)  :: con_xlv
        real(kind=kind_phys),intent(in)  :: con_xlf
        real(kind=kind_phys),intent(in)  :: con_rv
        real(kind=kind_phys),intent(in)  :: con_rd
        real(kind=kind_phys),intent(in)  :: con_ep2
        real(kind=kind_phys),intent(in)  :: con_grav
        real(kind=kind_phys),intent(in)  :: con_cp
        real(kind=kind_phys),intent(in)  :: con_cpv
        real(kind=kind_phys),intent(in)  :: con_rcp
        real(kind=kind_phys),intent(in)  :: con_p608
        real(kind=kind_phys),intent(in)  :: con_cliq
        real(kind=kind_phys),intent(in)  :: con_cice
        real(kind=kind_phys),intent(in)  :: con_karman
        real(kind=kind_phys),intent(in)  :: con_t0c

        ! Initialize CCPP error handling variables
        errmsg = ''
        errflg = 0

        xlv    = con_xlv
        xlf    = con_xlf
        r_v    = con_rv
        r_d    = con_rd
        ep_2   = con_ep2
        grav   = con_grav
        cp     = con_cp
        cpv    = con_cpv
        rcp    = con_rcp
        p608   = con_p608
        cliq   = con_cliq
        cice   = con_cice
        karman = con_karman
        t0c    = con_t0c
       
        xls    = xlv+xlf      != 2.85E6 (J/kg) sublimation                                      
        rvovrd = r_v/r_d      != 1.608
        ep_3   = 1.-ep_2      != 0.378
        gtr    = grav/tref
        rk     = cp/r_d
        tv0    = p608*tref
        tv1    = (1.+p608)*tref
        xlscp  = (xlv+xlf)/cp
        xlvcp  = xlv/cp
        g_inv  = 1./grav

        ! Consistency checks
        if (.not. do_mynnedmf) then
          errmsg = 'Logic error: do_mynnedmf = .false.'
          errflg = 1
            return
         end if

        if (lheatstrg) then
          errmsg = 'Logic error: lheatstrg not implemented for MYNN PBL'
          errflg = 1
          return
        end if

      end subroutine mynnedmf_wrapper_init

!>\defgroup gp_mynnedmf MYNN-EDMF PBL and Shallow Convection Module  
!> This scheme (1) performs pre-mynnedmf work, (2) runs the mynnedmf, and (3) performs post-mynnedmf work
!> \section arg_table_mynnedmf_wrapper_run Argument Table
!! \htmlinclude mynnedmf_wrapper_run.html
!!
SUBROUTINE mynnedmf_wrapper_run(        &
     &  im,levs,                        &
     &  flag_init,flag_restart,         &
     &  lssav, ldiag3d, qdiag3d,        &
     &  lsidea, cplflx,                 &
     &  delt,dtf,dx,zorl,               &
     &  phii,u,v,omega,t3d,             &
     &  qgrs_water_vapor,               &
     &  qgrs_liquid_cloud,              &
     &  qgrs_ice_cloud,                 &
     &  qgrs_cloud_droplet_num_conc,    &
     &  qgrs_cloud_ice_num_conc,        &
     &  qgrs_ozone,                     &
     &  qgrs_water_aer_num_conc,        &
     &  qgrs_ice_aer_num_conc,          &
     &  qgrs_cccn,                      &
     &  prsl,exner,                     &
     &  slmsk,tsurf,qsfc,ps,            &
     &  ust,ch,hflx,qflx,wspd,rb,       &
     &  dtsfc1,dqsfc1,                  &
     &  dusfc1,dvsfc1,                  &
     &  dusfci_diag,dvsfci_diag,        &
     &  dtsfci_diag,dqsfci_diag,        &
     &  dusfc_diag,dvsfc_diag,          &
     &  dtsfc_diag,dqsfc_diag,          &
     &  dusfc_cice,dvsfc_cice,          &
     &  dtsfc_cice,dqsfc_cice,          &
     &  hflx_wat,qflx_wat,stress_wat,   &
     &  oceanfrac,fice,wet,icy,dry,     &
     &  dusfci_cpl,dvsfci_cpl,          &
     &  dtsfci_cpl,dqsfci_cpl,          &
     &  dusfc_cpl,dvsfc_cpl,            &
     &  dtsfc_cpl,dqsfc_cpl,            &
     &  recmol,                         &
     &  qke,qke_adv,Tsq,Qsq,Cov,        &
     &  el_pbl,sh3d,sm3d,exch_h,exch_m, &
     &  dqke,qwt,qshear,qbuoy,qdiss,    &
     &  Pblh,kpbl,                      &
     &  qc_bl,qi_bl,cldfra_bl,          &
     &  edmf_a,edmf_w,edmf_qt,          &
     &  edmf_thl,edmf_ent,edmf_qc,      &
     &  sub_thl,sub_sqv,det_thl,det_sqv,&
     &  nupdraft,maxMF,ktop_plume,      &
     &  dudt, dvdt, dtdt,                                  &
     &  dqdt_water_vapor,            dqdt_liquid_cloud,    & ! <=== ntqv, ntcw
     &  dqdt_ice_cloud,              dqdt_ozone,           & ! <=== ntiw, ntoz
     &  dqdt_cloud_droplet_num_conc, dqdt_ice_num_conc,    & ! <=== ntlnc, ntinc
     &  dqdt_water_aer_num_conc,     dqdt_ice_aer_num_conc,& ! <=== ntwa, ntia
     &  dqdt_cccn,                                         & ! <=== ntccn
     &  flag_for_pbl_generic_tend,                         &
     &  dtend, dtidx, index_of_temperature,                &
     &  index_of_x_wind, index_of_y_wind, ntke,            &
     &  ntqv, ntcw, ntiw, ntoz, ntlnc, ntinc, ntwa, ntia,  &
     &  index_of_process_pbl, htrsw, htrlw, xmu,           &
     &  bl_mynn_tkebudget,     bl_mynn_tkeadvect,          &
     &  bl_mynn_cloudpdf,      bl_mynn_mixlength,          &
     &  bl_mynn_edmf,                                      &
     &  bl_mynn_edmf_mom,      bl_mynn_edmf_tke,           &
     &  bl_mynn_cloudmix,      bl_mynn_mixqt,              &
     &  bl_mynn_output,        bl_mynn_closure,            &
     &  icloud_bl, do_mynnsfclay,                          &
     &  imp_physics, imp_physics_gfdl,                     &
     &  imp_physics_thompson, imp_physics_wsm6,            &
     &  chem3d, frp, mix_chem, rrfs_smoke, fire_turb, nchem, ndvel, &
     &  imp_physics_nssl, nssl_ccn_on,                     &
     &  ltaerosol, mraerosol, spp_wts_pbl, spp_pbl, lprnt, huge, errmsg, errflg  )

! should be moved to inside the mynn:
     use machine,        only: kind_phys
     use bl_mynn_common, only: cp, r_d, grav, g_inv, zero, &
         xlv, xlvcp, xlscp
     use module_bl_mynn, only: mynn_bl_driver

!------------------------------------------------------------------- 
     implicit none
!------------------------------------------------------------------- 

     real(kind=kind_phys)          :: huge
     character(len=*), intent(out) :: errmsg
     integer, intent(out)          :: errflg

     logical, intent(in) :: lssav, ldiag3d, lsidea, qdiag3d
     logical, intent(in) :: cplflx

     !smoke/chem
     integer, intent(in) :: nchem, ndvel
     integer, parameter  :: kdvel=1

! NAMELIST OPTIONS (INPUT):
     logical, intent(in) ::                                 &
     &       bl_mynn_tkeadvect,                             &
     &       bl_mynn_tkebudget,                             &
     &       ltaerosol, mraerosol,                          &
     &       lprnt,                                         &
     &       do_mynnsfclay,                                 &
     &       flag_for_pbl_generic_tend,                     &
     &       nssl_ccn_on
      integer, intent(in) ::                                &
     &       bl_mynn_cloudpdf,                              &
     &       bl_mynn_mixlength,                             &
     &       icloud_bl,                                     &
     &       bl_mynn_edmf,                                  &
     &       bl_mynn_edmf_mom,                              &
     &       bl_mynn_edmf_tke,                              &
     &       bl_mynn_cloudmix,                              &
     &       bl_mynn_mixqt,                                 &
     &       bl_mynn_output,                                &
     &       imp_physics, imp_physics_wsm6,                 &
     &       imp_physics_thompson, imp_physics_gfdl,        &
     &       imp_physics_nssl,                              &
     &       spp_pbl
      real, intent(in) ::                                   &
     &       bl_mynn_closure

!TENDENCY DIAGNOSTICS
      real(kind=kind_phys), intent(inout), optional :: dtend(:,:,:)
      integer, intent(in) :: dtidx(:,:)
      integer, intent(in) :: index_of_temperature, index_of_x_wind
      integer, intent(in) :: index_of_y_wind, index_of_process_pbl
      integer, intent(in) :: ntoz, ntqv, ntcw, ntiw, ntlnc
      integer, intent(in) :: ntinc, ntwa, ntia, ntke

!MISC CONFIGURATION OPTIONS
      INTEGER, PARAMETER ::                                 &
     &       bl_mynn_mixscalars=1
      LOGICAL ::                                            &
     &       FLAG_QI, FLAG_QNI, FLAG_QC, FLAG_QNC,          &
     &       FLAG_QNWFA, FLAG_QNIFA, FLAG_OZONE
      ! Define locally until needed from CCPP
      LOGICAL, PARAMETER :: cycling = .false.
      INTEGER, PARAMETER :: param_first_scalar = 1
      INTEGER ::                                            &
     &      p_qc, p_qr, p_qi, p_qs, p_qg, p_qnc, p_qni

!MYNN-1D
      REAL(kind=kind_phys), intent(in) :: delt, dtf
      INTEGER, intent(in) :: im, levs
      LOGICAL, intent(in) :: flag_init, flag_restart
      INTEGER :: initflag, k, i
      INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE,                   &
     &           IMS,IME,JMS,JME,KMS,KME,                   &
     &           ITS,ITE,JTS,JTE,KTS,KTE

      REAL(kind=kind_phys) :: tem

!MYNN-3D
      real(kind=kind_phys), dimension(:,:), intent(in)    :: phii
      real(kind=kind_phys), dimension(:,:), intent(inout) ::             &
     &        dtdt, dudt, dvdt,                                          &
     &        dqdt_water_vapor, dqdt_liquid_cloud, dqdt_ice_cloud,       &
     &        dqdt_cloud_droplet_num_conc, dqdt_ice_num_conc,            &
     &        dqdt_ozone, dqdt_water_aer_num_conc, dqdt_ice_aer_num_conc
      real(kind=kind_phys), dimension(:,:), intent(inout) ::dqdt_cccn
      real(kind=kind_phys), dimension(:,:), intent(inout) ::             &
     &        qke, qke_adv, EL_PBL, Sh3D, Sm3D,                          &
     &        qc_bl, qi_bl, cldfra_bl
     !These 10 arrays are only allocated when bl_mynn_output > 0
      real(kind=kind_phys), dimension(:,:), intent(inout) ::             &
     &        edmf_a,edmf_w,edmf_qt,                                     &
     &        edmf_thl,edmf_ent,edmf_qc,                                 &
     &        sub_thl,sub_sqv,det_thl,det_sqv
      real(kind=kind_phys), dimension(:,:), intent(inout) ::             &
     &        dqke,qWT,qSHEAR,qBUOY,qDISS
      real(kind=kind_phys), dimension(:,:), intent(inout) ::             &
     &        t3d,qgrs_water_vapor,qgrs_liquid_cloud,qgrs_ice_cloud
      real(kind=kind_phys), dimension(:,:), intent(in) ::                &
     &        u,v,omega,                                                 &
     &        exner,prsl,                                                &
     &        qgrs_cloud_droplet_num_conc,                               &
     &        qgrs_cloud_ice_num_conc,                                   &
     &        qgrs_ozone,                                                &
     &        qgrs_water_aer_num_conc,                                   &
     &        qgrs_ice_aer_num_conc
      real(kind=kind_phys), dimension(:,:), intent(in) ::qgrs_cccn
      real(kind=kind_phys), dimension(:,:), intent(out) ::               &
     &        Tsq, Qsq, Cov, exch_h, exch_m
      real(kind=kind_phys), dimension(:), intent(in) :: xmu
      real(kind=kind_phys), dimension(:,:), intent(in) :: htrsw, htrlw
      ! spp_wts_pbl only allocated if spp_pbl == 1
      real(kind_phys), dimension(:,:),       intent(in) :: spp_wts_pbl

     !LOCAL
      real(kind=kind_phys), dimension(im,levs) ::                        &
     &        sqv,sqc,sqi,qnc,qni,ozone,qnwfa,qnifa,                     &
     &        dz, w, p, rho, th, qv, delp,                               &
     &        RUBLTEN, RVBLTEN, RTHBLTEN, RQVBLTEN,                      &
     &        RQCBLTEN, RQNCBLTEN, RQIBLTEN, RQNIBLTEN,                  &
     &        RQNWFABLTEN, RQNIFABLTEN
      real(kind=kind_phys), allocatable :: old_ozone(:,:)

!smoke/chem arrays
      real(kind_phys), dimension(:), intent(inout) :: frp
      logical, intent(in) :: mix_chem, fire_turb, rrfs_smoke
      real(kind=kind_phys), dimension(:,:,:), intent(inout) :: chem3d
      real(kind=kind_phys), dimension(im)   :: emis_ant_no
      real(kind=kind_phys), dimension(im,ndvel) :: vdep

!MYNN-2D
      real(kind=kind_phys), dimension(:), intent(in) ::                  &
     &        dx,zorl,slmsk,tsurf,qsfc,ps,                               &
     &        hflx,qflx,ust,wspd,rb,recmol

      real(kind=kind_phys), dimension(:), intent(in) ::                  &
     &        dusfc_cice,dvsfc_cice,dtsfc_cice,dqsfc_cice,               &
     &        stress_wat,hflx_wat,qflx_wat,                              &
     &        oceanfrac,fice

      logical, dimension(:), intent(in) ::                               &
     &        wet, dry, icy

      real(kind=kind_phys), dimension(:), intent(inout) ::               &
     &        pblh,dusfc_diag,dvsfc_diag,dtsfc_diag,dqsfc_diag
      real(kind=kind_phys), dimension(:), intent(out) ::                 &
     &        ch,dtsfc1,dqsfc1,dusfc1,dvsfc1,                            &
     &        dtsfci_diag,dqsfci_diag,dusfci_diag,dvsfci_diag,           &
     &        maxMF
      integer, dimension(:), intent(inout) ::                            &
     &        kpbl,nupdraft,ktop_plume

      real(kind=kind_phys), dimension(:), intent(inout) ::               &
     &        dusfc_cpl,dvsfc_cpl,dtsfc_cpl,dqsfc_cpl
      real(kind=kind_phys), dimension(:), intent(inout) ::               &
     &        dusfci_cpl,dvsfci_cpl,dtsfci_cpl,dqsfci_cpl

     !LOCAL
      real, dimension(im) ::                                             &
     &        hfx,qfx,rmol,xland,uoce,voce,vdfg,znt,ts
      integer :: idtend
      real, dimension(im) :: dusfci1,dvsfci1,dtsfci1,dqsfci1
      real(kind=kind_phys), allocatable :: save_qke_adv(:,:)

      ! Initialize CCPP error handling variables
      errmsg = ''
      errflg = 0

      if (lprnt) then
         write(0,*)"=============================================="
         write(0,*)"in mynn wrapper..."
         write(0,*)"flag_init=",flag_init
         write(0,*)"flag_restart=",flag_restart
      endif

      if (.not. flag_for_pbl_generic_tend .and. ldiag3d) then
         idtend = dtidx(ntke+100,index_of_process_pbl)
         if (idtend>=1) then
            allocate(save_qke_adv(im,levs))
            save_qke_adv=qke_adv
         endif
      endif

      ! DH* TODO: Use flag_restart to distinguish which fields need
      ! to be initialized and which are read from restart files
      if (flag_init) then
         initflag=1
         !print*,"in MYNN, initflag=",initflag
      else
         initflag=0
         !print*,"in MYNN, initflag=",initflag
      endif

      !initialize arrays for test
      EMIS_ANT_NO = 0.
      vdep = 0. ! hli for chem dry deposition, 0 temporarily

  ! Check incoming moist species to ensure non-negative values
  ! First, create height (dz) and pressure differences (delp) 
  ! across model layers
      do k=1,levs
         do i=1,im
            dz(i,k)=(phii(i,k+1) - phii(i,k))*g_inv
         enddo
      enddo

      do i=1,im
         delp(i,1)  = ps(i) - (prsl(i,2)*dz(i,1) + prsl(i,1)*dz(i,2))/(dz(i,1)+dz(i,2))
         do k=2,levs-1
            delp(i,k) = (prsl(i,k)*dz(i,k-1) + prsl(i,k-1)*dz(i,k))/(dz(i,k)+dz(i,k-1)) - &
                        (prsl(i,k+1)*dz(i,k) + prsl(i,k)*dz(i,k+1))/(dz(i,k)+dz(i,k+1))
         enddo
         delp(i,levs) = delp(i,levs-1)
      enddo

      do i=1,im
         call moisture_check2(levs, delt,            &
                              delp(i,:), exner(i,:), &
                              qgrs_water_vapor(i,:), &
                              qgrs_liquid_cloud(i,:),&
                              qgrs_ice_cloud(i,:),   &
                              t3d(i,:)               )
      enddo

      FLAG_OZONE = ntoz>0

  ! Assign variables for each microphysics scheme
        if (imp_physics == imp_physics_wsm6) then
  ! WSM6
         FLAG_QI = .true.
         FLAG_QNI= .false.
         FLAG_QC = .true.
         FLAG_QNC= .false.
         FLAG_QNWFA= .false.
         FLAG_QNIFA= .false.
         p_qc = 2
         p_qr = 0
         p_qi = 2 
         p_qs = 0 
         p_qg = 0
         p_qnc= 0 
         p_qni= 0 
         do k=1,levs
            do i=1,im
              sqv(i,k)   = qgrs_water_vapor(i,k)
              sqc(i,k)   = qgrs_liquid_cloud(i,k)
              sqi(i,k)   = qgrs_ice_cloud(i,k)
              ozone(i,k) = qgrs_ozone(i,k)
              qnc(i,k)   = 0.
              qni(i,k)   = 0.
              qnwfa(i,k) = 0.
              qnifa(i,k) = 0.
            enddo
          enddo
        elseif (imp_physics == imp_physics_nssl ) then
  ! NSSL
         FLAG_QI = .true.
         FLAG_QNI= .true.
         FLAG_QC = .true.
         FLAG_QNC= .true.
         FLAG_QNWFA= nssl_ccn_on ! ERM: Perhaps could use this field for CCN field?
         FLAG_QNIFA= .false.
         ! p_q vars not used?
         p_qc = 2
         p_qr = 0
         p_qi = 2 
         p_qs = 0 
         p_qg = 0
         p_qnc= 0 
         p_qni= 0 
         do k=1,levs
            do i=1,im
              sqv(i,k)  = qgrs_water_vapor(i,k)
              sqc(i,k)    = qgrs_liquid_cloud(i,k)
              sqi(i,k)    = qgrs_ice_cloud(i,k)
              ozone(i,k) = qgrs_ozone(i,k)
              qnc(i,k)   = qgrs_cloud_droplet_num_conc(i,k)
              qni(i,k)   = qgrs_cloud_ice_num_conc(i,k)
              qnwfa(i,k) = 0.
              IF ( nssl_ccn_on ) THEN
                qnwfa(i,k) = qgrs_cccn(i,k)
              ENDIF
              qnifa(i,k) = 0.
            enddo
          enddo
        elseif (imp_physics == imp_physics_thompson) then
  ! Thompson
          if(ltaerosol) then
            FLAG_QI = .true.
            FLAG_QNI= .true.
            FLAG_QC = .true.
            FLAG_QNC= .true.
            FLAG_QNWFA= .true.
            FLAG_QNIFA= .true.
            p_qc = 2
            p_qr = 0
            p_qi = 2
            p_qs = 0
            p_qg = 0
            p_qnc= 0
            p_qni= 0
            do k=1,levs
              do i=1,im
                sqv(i,k)   = qgrs_water_vapor(i,k)
                sqc(i,k)   = qgrs_liquid_cloud(i,k)
                sqi(i,k)   = qgrs_ice_cloud(i,k)
                qnc(i,k)   = qgrs_cloud_droplet_num_conc(i,k)
                qni(i,k)   = qgrs_cloud_ice_num_conc(i,k)
                ozone(i,k) = qgrs_ozone(i,k)
                qnwfa(i,k) = qgrs_water_aer_num_conc(i,k)
                qnifa(i,k) = qgrs_ice_aer_num_conc(i,k)
              enddo
            enddo
          else if(mraerosol) then
            FLAG_QI = .true.
            FLAG_QNI= .true.
            FLAG_QC = .true.
            FLAG_QNC= .true.
            FLAG_QNWFA= .false.
            FLAG_QNIFA= .false.
            p_qc = 2
            p_qr = 0
            p_qi = 2
            p_qs = 0
            p_qg = 0
            p_qnc= 0
            p_qni= 0
            do k=1,levs
              do i=1,im
                sqv(i,k)   = qgrs_water_vapor(i,k)
                sqc(i,k)   = qgrs_liquid_cloud(i,k)
                sqi(i,k)   = qgrs_ice_cloud(i,k)
                qnc(i,k)   = qgrs_cloud_droplet_num_conc(i,k)
                qni(i,k)   = qgrs_cloud_ice_num_conc(i,k)
                ozone(i,k) = qgrs_ozone(i,k)
                qnwfa(i,k) = 0.
                qnifa(i,k) = 0.
              enddo
            enddo
          else
            FLAG_QI = .true.
            FLAG_QNI= .true.
            FLAG_QC = .true.
            FLAG_QNC= .false.
            FLAG_QNWFA= .false.
            FLAG_QNIFA= .false.
            p_qc = 2
            p_qr = 0
            p_qi = 2
            p_qs = 0
            p_qg = 0
            p_qnc= 0
            p_qni= 0
            do k=1,levs
              do i=1,im
                sqv(i,k)   = qgrs_water_vapor(i,k)
                sqc(i,k)   = qgrs_liquid_cloud(i,k)
                sqi(i,k)   = qgrs_ice_cloud(i,k)
                qnc(i,k)   = 0.
                qni(i,k)   = qgrs_cloud_ice_num_conc(i,k)
                ozone(i,k) = qgrs_ozone(i,k)
                qnwfa(i,k) = 0.
                qnifa(i,k) = 0.
              enddo
            enddo
          endif
        elseif (imp_physics == imp_physics_gfdl) then
  ! GFDL MP
          FLAG_QI = .true.
          FLAG_QNI= .false.
          FLAG_QC = .true.
          FLAG_QNC= .false.
          FLAG_QNWFA= .false.
          FLAG_QNIFA= .false.
          p_qc = 2
          p_qr = 0
          p_qi = 2
          p_qs = 0
          p_qg = 0
          p_qnc= 0
          p_qni= 0
          do k=1,levs
            do i=1,im
                sqv(i,k)   = qgrs_water_vapor(i,k)
                sqc(i,k)   = qgrs_liquid_cloud(i,k)
                sqi(i,k)   = qgrs_ice_cloud(i,k)
                qnc(i,k)   = 0.
                qni(i,k)   = 0.
                qnwfa(i,k) = 0.
                qnifa(i,k) = 0.
                ozone(i,k) = qgrs_ozone(i,k)
            enddo
          enddo
        else
          print*,"In MYNN wrapper. Unknown microphysics scheme, imp_physics=",imp_physics
          print*,"Defaulting to qc and qv species only..."
          FLAG_QI = .false.
          FLAG_QNI= .false.
          FLAG_QC = .true.
          FLAG_QNC= .false.
          FLAG_QNWFA= .false.
          FLAG_QNIFA= .false.
          p_qc = 2
          p_qr = 0
          p_qi = 0
          p_qs = 0
          p_qg = 0
          p_qnc= 0
          p_qni= 0
          do k=1,levs
            do i=1,im
                sqv(i,k)   = qgrs_water_vapor(i,k)
                sqc(i,k)   = qgrs_liquid_cloud(i,k)
                sqi(i,k)   = 0.
                qnc(i,k)   = 0.
                qni(i,k)   = 0.
                qnwfa(i,k) = 0.
                qnifa(i,k) = 0.
                ozone(i,k) = qgrs_ozone(i,k)
            enddo
          enddo
        endif
       if(ldiag3d .and. dtidx(100+ntoz,index_of_process_pbl)>1) then
         allocate(old_ozone(im,levs))
         old_ozone = ozone
       endif
       if (lprnt)write(0,*)"prepping MYNN-EDMF variables..."

       do k=1,levs
          do i=1,im
          !   dz(i,k)=(phii(i,k+1) - phii(i,k))*g_inv
             th(i,k)=t3d(i,k)/exner(i,k)
          ! keep as specific humidity
          !   qv(i,k)=qvsh(i,k)/(1.0 - qvsh(i,k))
          !   qc(i,k)=qc(i,k)/(1.0 - qvsh(i,k))
          !   qi(i,k)=qi(i,k)/(1.0 - qvsh(i,k))
             rho(i,k)=prsl(i,k)/(r_d*t3d(i,k))
             w(i,k) = -omega(i,k)/(rho(i,k)*grav)
         enddo
      enddo

      do i=1,im
         if (slmsk(i)==1. .or. slmsk(i)==2.) then !sea/land/ice mask (=0/1/2) in FV3
            xland(i)=1.0                          !but land/water = (1/2) in SFCLAY_mynn
         else
            xland(i)=2.0
         endif
         uoce(i)=0.0
         voce(i)=0.0
         vdfg(i)=0.0
         !ust(i) = sqrt(stress(i))
         ch(i)=0.0
         hfx(i)=hflx(i)*rho(i,1)*cp
         qfx(i)=qflx(i)*rho(i,1)

         dtsfc1(i) = hfx(i)
         dqsfc1(i) = qfx(i)*XLV
         dusfc1(i) = -1.*rho(i,1)*ust(i)*ust(i)*u(i,1)/wspd(i)
         dvsfc1(i) = -1.*rho(i,1)*ust(i)*ust(i)*v(i,1)/wspd(i)

         !BWG: diagnostic surface fluxes for scalars & momentum
         dtsfci_diag(i) = dtsfc1(i)
         dqsfci_diag(i) = dqsfc1(i)
         dtsfc_diag(i)  = dtsfc_diag(i) + dtsfc1(i)*delt
         dqsfc_diag(i)  = dqsfc_diag(i) + dqsfc1(i)*delt
         dusfci_diag(i) = dusfc1(i)
         dvsfci_diag(i) = dvsfc1(i)
         dusfc_diag(i)  = dusfc_diag(i) + dusfci_diag(i)*delt
         dvsfc_diag(i)  = dvsfc_diag(i) + dvsfci_diag(i)*delt

         znt(i)=zorl(i)*0.01 !cm -> m?
         if (do_mynnsfclay) then
           rmol(i)=recmol(i)
         else
           if (hfx(i) .ge. 0.)then
             rmol(i)=-hfx(i)/(200.*dz(i,1)*0.5)
           else
             rmol(i)=ABS(rb(i))*1./(dz(i,1)*0.5)
           endif
         endif
         ts(i)=tsurf(i)/exner(i,1)  !theta
!        qsfc(i)=qss(i)
!        ps(i)=pgr(i)
!        wspd(i)=wind(i)
      enddo

      ! BWG: Coupling insertion
      if (cplflx) then
        do i=1,im
          if (oceanfrac(i) > zero) then ! Ocean only, NO LAKES
            if ( .not. wet(i)) then ! no open water, use results from CICE
              dusfci_cpl(i) = dusfc_cice(i)
              dvsfci_cpl(i) = dvsfc_cice(i)
              dtsfci_cpl(i) = dtsfc_cice(i)
              dqsfci_cpl(i) = dqsfc_cice(i)
            elseif (icy(i) .or. dry(i)) then ! use stress_ocean for opw component at mixed point
              if (wspd(i) > zero) then
                dusfci_cpl(i) = -1.*rho(i,1)*stress_wat(i)*u(i,1)/wspd(i)   ! U-momentum flux
                dvsfci_cpl(i) = -1.*rho(i,1)*stress_wat(i)*v(i,1)/wspd(i)   ! V-momentum flux
              else
                dusfci_cpl(i) = zero
                dvsfci_cpl(i) = zero
              endif
              dtsfci_cpl(i) =  cp*rho(i,1)*hflx_wat(i) ! sensible heat flux over open ocean
              dqsfci_cpl(i) = XLV*rho(i,1)*qflx_wat(i) ! latent heat flux over open ocean
            else                                       ! use results from this scheme for 100% open ocean
              dusfci_cpl(i) = dusfci_diag(i)
              dvsfci_cpl(i) = dvsfci_diag(i)
              dtsfci_cpl(i) = dtsfci_diag(i)
              dqsfci_cpl(i) = dqsfci_diag(i)
            endif
!
            dusfc_cpl (i) = dusfc_cpl(i) + dusfci_cpl(i) * delt
            dvsfc_cpl (i) = dvsfc_cpl(i) + dvsfci_cpl(i) * delt
            dtsfc_cpl (i) = dtsfc_cpl(i) + dtsfci_cpl(i) * delt
            dqsfc_cpl (i) = dqsfc_cpl(i) + dqsfci_cpl(i) * delt
          else ! If no ocean
            dusfc_cpl(i) = huge
            dvsfc_cpl(i) = huge
            dtsfc_cpl(i) = huge
            dqsfc_cpl(i) = huge
          endif ! Ocean only, NO LAKES
        enddo
      endif ! End coupling insertion

      if (lprnt) then
         print*
         write(0,*)"===CALLING mynn_bl_driver; input:"
         print*,"bl_mynn_tkebudget=",bl_mynn_tkebudget," bl_mynn_tkeadvect=",bl_mynn_tkeadvect
         print*,"bl_mynn_cloudpdf=",bl_mynn_cloudpdf," bl_mynn_mixlength=",bl_mynn_mixlength
         print*,"bl_mynn_edmf=",bl_mynn_edmf," bl_mynn_edmf_mom=",bl_mynn_edmf_mom
         print*,"bl_mynn_edmf_tke=",bl_mynn_edmf_tke
         print*,"bl_mynn_cloudmix=",bl_mynn_cloudmix," bl_mynn_mixqt=",bl_mynn_mixqt
         print*,"icloud_bl=",icloud_bl
         print*,"T:",t3d(1,1),t3d(1,2),t3d(1,levs)
         print*,"TH:",th(1,1),th(1,2),th(1,levs)
         print*,"rho:",rho(1,1),rho(1,2),rho(1,levs)
         print*,"exner:",exner(1,1),exner(1,2),exner(1,levs)
         print*,"prsl:",prsl(1,1),prsl(1,2),prsl(1,levs)
         print*,"dz:",dz(1,1),dz(1,2),dz(1,levs)
         print*,"u:",u(1,1),u(1,2),u(1,levs)
         print*,"v:",v(1,1),v(1,2),v(1,levs)
         print*,"sqv:",sqv(1,1),sqv(1,2),sqv(1,levs)
         print*,"sqc:",sqc(1,1),sqc(1,2),sqc(1,levs)
         print*,"sqi:",sqi(1,1),sqi(1,2),sqi(1,levs)
         print*,"rmol:",rmol(1)," ust:",ust(1)
         print*," dx=",dx(1),"initflag=",initflag
         print*,"Tsurf:",tsurf(1)," Thetasurf:",ts(1)
         print*,"HFX:",hfx(1)," qfx",qfx(1)
         print*,"qsfc:",qsfc(1)," ps:",ps(1)
         print*,"wspd:",wspd(1)," rb=",rb(1)
         print*,"znt:",znt(1)," delt=",delt
         print*,"im=",im," levs=",levs
         print*,"PBLH=",pblh(1)," KPBL=",KPBL(1)," xland=",xland(1)
         print*,"vdfg=",vdfg(1)," ch=",ch(1)
         !print*,"TKE:",TKE_PBL(1,1),TKE_PBL(1,2),TKE_PBL(1,levs)
         print*,"qke:",qke(1,1),qke(1,2),qke(1,levs)
         print*,"el_pbl:",el_pbl(1,1),el_pbl(1,2),el_pbl(1,levs)
         print*,"Sh3d:",Sh3d(1,1),sh3d(1,2),sh3d(1,levs)
         !print*,"exch_h:",exch_h(1,1),exch_h(1,2),exch_h(1,levs) ! - intent(out)
         !print*,"exch_m:",exch_m(1,1),exch_m(1,2),exch_m(1,levs) ! - intent(out)
         print*,"max cf_bl:",maxval(cldfra_bl(1,:))
      endif


              CALL  mynn_bl_driver(                                    &
     &             initflag=initflag,restart=flag_restart,             &
     &             cycling=cycling,                                    &
     &             delt=delt,dz=dz,dx=dx,znt=znt,                      &
     &             u=u,v=v,w=w,th=th,sqv3D=sqv,sqc3D=sqc,              &
     &             sqi3D=sqi,qnc=qnc,qni=qni,                          &
     &             qnwfa=qnwfa,qnifa=qnifa,ozone=ozone,                &
     &             p=prsl,exner=exner,rho=rho,T3D=t3d,                 &
     &             xland=xland,ts=ts,qsfc=qsfc,ps=ps,                  &
     &             ust=ust,ch=ch,hfx=hfx,qfx=qfx,rmol=rmol,            &
     &             wspd=wspd,uoce=uoce,voce=voce,vdfg=vdfg,            & !input
     &             qke=QKE,qke_adv=qke_adv,                            & !output
     &             sh3d=Sh3d,sm3d=Sm3d,                                &
!chem/smoke
     &             nchem=nchem,kdvel=kdvel,ndvel=ndvel,                &
     &             Chem3d=chem3d,Vdep=vdep,                            &
     &             FRP=frp,EMIS_ANT_NO=emis_ant_no,                    &
     &             mix_chem=mix_chem,fire_turb=fire_turb,              &
     &             rrfs_smoke=rrfs_smoke,                              &
!-----
     &             Tsq=tsq,Qsq=qsq,Cov=cov,                            & !output
     &             RUBLTEN=RUBLTEN,RVBLTEN=RVBLTEN,RTHBLTEN=RTHBLTEN,  & !output
     &             RQVBLTEN=RQVBLTEN,RQCBLTEN=rqcblten,                &
     &             RQIBLTEN=rqiblten,RQNCBLTEN=rqncblten,              & !output
     &             RQNIBLTEN=rqniblten,RQNWFABLTEN=RQNWFABLTEN,        & !output
     &             RQNIFABLTEN=RQNIFABLTEN,dozone=dqdt_ozone,          & !output
     &             EXCH_H=exch_h,EXCH_M=exch_m,                        & !output
     &             pblh=pblh,KPBL=KPBL,                                & !output
     &             el_pbl=el_pbl,                                      & !output
     &             dqke=dqke,                                          & !output
     &             qWT=qWT,qSHEAR=qSHEAR,qBUOY=qBUOY,qDISS=qDISS,      & !output
     &             bl_mynn_tkeadvect=bl_mynn_tkeadvect,                &
     &             bl_mynn_tkebudget=bl_mynn_tkebudget,                & !input parameter
     &             bl_mynn_cloudpdf=bl_mynn_cloudpdf,                  & !input parameter
     &             bl_mynn_mixlength=bl_mynn_mixlength,                & !input parameter
     &             icloud_bl=icloud_bl,                                & !input parameter
     &             qc_bl=qc_bl,qi_bl=qi_bl,cldfra_bl=cldfra_bl,        & !output
     &             closure=bl_mynn_closure,bl_mynn_edmf=bl_mynn_edmf,  & !input parameter
     &             bl_mynn_edmf_mom=bl_mynn_edmf_mom,                  & !input parameter
     &             bl_mynn_edmf_tke=bl_mynn_edmf_tke,                  & !input parameter
     &             bl_mynn_mixscalars=bl_mynn_mixscalars,              & !input parameter
     &             bl_mynn_output=bl_mynn_output,                      & !input parameter
     &             bl_mynn_cloudmix=bl_mynn_cloudmix,                  & !input parameter
     &             bl_mynn_mixqt=bl_mynn_mixqt,                        & !input parameter
     &             edmf_a=edmf_a,edmf_w=edmf_w,edmf_qt=edmf_qt,        & !output
     &             edmf_thl=edmf_thl,edmf_ent=edmf_ent,edmf_qc=edmf_qc,&!output
     &             sub_thl3D=sub_thl,sub_sqv3D=sub_sqv,                &
     &             det_thl3D=det_thl,det_sqv3D=det_sqv,                &
     &             nupdraft=nupdraft,maxMF=maxMF,                      & !output
     &             ktop_plume=ktop_plume,                              & !output
     &             spp_pbl=spp_pbl,pattern_spp_pbl=spp_wts_pbl,        & !input
     &             RTHRATEN=htrlw,                                     & !input
     &             FLAG_QI=flag_qi,FLAG_QNI=flag_qni,                  & !input
     &             FLAG_QC=flag_qc,FLAG_QNC=flag_qnc,                  & !input
     &             FLAG_QNWFA=FLAG_QNWFA,FLAG_QNIFA=FLAG_QNIFA,        & !input
     &             FLAG_OZONE=FLAG_OZONE,                              & !input
     &             IDS=1,IDE=im,JDS=1,JDE=1,KDS=1,KDE=levs,            & !input
     &             IMS=1,IME=im,JMS=1,JME=1,KMS=1,KME=levs,            & !input
     &             ITS=1,ITE=im,JTS=1,JTE=1,KTS=1,KTE=levs)              !input


     ! POST MYNN (INTERSTITIAL) WORK:
        !update/save MYNN-only variables
        !do k=1,levs
        !   do i=1,im
        !      gq0(i,k,4)=qke(i,k,1)      !tke*2
        !   enddo
        !enddo
        !For MYNN, convert TH-tend to T-tend
        do k = 1, levs
           do i = 1, im
              dtdt(i,k) = dtdt(i,k) + RTHBLTEN(i,k)*exner(i,k)
              dudt(i,k) = dudt(i,k) + RUBLTEN(i,k)
              dvdt(i,k) = dvdt(i,k) + RVBLTEN(i,k)
           enddo
        enddo
        accum_duvt3dt: if(ldiag3d .or. lsidea) then
          call dtend_helper(index_of_x_wind,RUBLTEN)
          call dtend_helper(index_of_y_wind,RVBLTEN)
          call dtend_helper(index_of_temperature,RTHBLTEN,exner)
          if(ldiag3d) then
            call dtend_helper(100+ntoz,dqdt_ozone)
            ! idtend = dtidx(100+ntoz,index_of_process_pbl)
            ! if(idtend>=1) then
            !   dtend(:,:,idtend) = dtend(:,:,idtend) + (ozone-old_ozone)
            !   deallocate(old_ozone)
            ! endif
          endif
        endif accum_duvt3dt
        !Update T, U and V:
        !do k = 1, levs
        !   do i = 1, im
        !      T3D(i,k) = T3D(i,k) + RTHBLTEN(i,k)*exner(i,k)*delt
        !      u(i,k)   = u(i,k) + RUBLTEN(i,k)*delt
        !      v(i,k)   = v(i,k) + RVBLTEN(i,k)*delt
        !   enddo
        !enddo

        !DO moist/scalar/tracer tendencies:
        if (imp_physics == imp_physics_wsm6) then
           ! WSM6
           do k=1,levs
             do i=1,im
               dqdt_water_vapor(i,k)  = RQVBLTEN(i,k) !/(1.0 + qv(i,k))
               dqdt_liquid_cloud(i,k) = RQCBLTEN(i,k) !/(1.0 + qv(i,k))
               dqdt_ice_cloud(i,k)    = RQIBLTEN(i,k) !/(1.0 + qv(i,k))
               !dqdt_ozone(i,k)        = 0.0
             enddo
           enddo
           if(ldiag3d .and. .not. flag_for_pbl_generic_tend) then
             call dtend_helper(100+ntqv,RQVBLTEN)
             call dtend_helper(100+ntcw,RQCBLTEN)
             call dtend_helper(100+ntiw,RQIBLTEN)
           endif
           !Update moist species:
           !do k=1,levs
           !  do i=1,im
           !    qgrs_water_vapor(i,k)  = qgrs_water_vapor(i,k)  + (RQVBLTEN(i,k)/(1.0+RQVBLTEN(i,k)))*delt
           !    qgrs_liquid_cloud(i,k) = qgrs_liquid_cloud(i,k) + RQCBLTEN(i,k)*delt
           !    qgrs_ice_cloud(i,k)    = qgrs_ice_cloud(i,k)    + RQIBLTEN(i,k)*delt
           !    !dqdt_ozone(i,k)        = 0.0
           !  enddo
           !enddo
        elseif (imp_physics == imp_physics_thompson) then
           ! Thompson-Aerosol
           if(ltaerosol) then
             do k=1,levs
               do i=1,im
                 dqdt_water_vapor(i,k)             = RQVBLTEN(i,k) !/(1.0 + qv(i,k))
                 dqdt_liquid_cloud(i,k)            = RQCBLTEN(i,k) !/(1.0 + qv(i,k))
                 dqdt_cloud_droplet_num_conc(i,k)  = RQNCBLTEN(i,k)
                 dqdt_ice_cloud(i,k)               = RQIBLTEN(i,k) !/(1.0 + qv(i,k))
                 dqdt_ice_num_conc(i,k)            = RQNIBLTEN(i,k)
                 !dqdt_ozone(i,k)                   = 0.0
                 dqdt_water_aer_num_conc(i,k)      = RQNWFABLTEN(i,k)
                 dqdt_ice_aer_num_conc(i,k)        = RQNIFABLTEN(i,k)
               enddo
             enddo
             if(ldiag3d .and. .not. flag_for_pbl_generic_tend) then
               call dtend_helper(100+ntqv,RQVBLTEN)
               call dtend_helper(100+ntcw,RQCBLTEN)
               call dtend_helper(100+ntlnc,RQNCBLTEN)
               call dtend_helper(100+ntiw,RQIBLTEN)
               call dtend_helper(100+ntinc,RQNIBLTEN)
               call dtend_helper(100+ntwa,RQNWFABLTEN)
               call dtend_helper(100+ntia,RQNIFABLTEN)
             endif
             !do k=1,levs
             !  do i=1,im
             !    qgrs_water_vapor(i,k)            = qgrs_water_vapor(i,k)    + (RQVBLTEN(i,k)/(1.0+RQVBLTEN(i,k)))*delt
             !    qgrs_liquid_cloud(i,k)           = qgrs_liquid_cloud(i,k)   + RQCBLTEN(i,k)*delt
             !    qgrs_ice_cloud(i,k)              = qgrs_ice_cloud(i,k)      + RQIBLTEN(i,k)*delt
             !    qgrs_cloud_droplet_num_conc(i,k) = qgrs_cloud_droplet_num_conc(i,k) + RQNCBLTEN(i,k)*delt
             !    qgrs_cloud_ice_num_conc(i,k)     = qgrs_cloud_ice_num_conc(i,k)     + RQNIBLTEN(i,k)*delt
             !    !dqdt_ozone(i,k)        = 0.0
             !    !qgrs_water_aer_num_conc(i,k)     = qgrs_water_aer_num_conc(i,k)     + RQNWFABLTEN(i,k)*delt
             !    !qgrs_ice_aer_num_conc(i,k)       = qgrs_ice_aer_num_conc(i,k)       + RQNIFABLTEN(i,k)*delt
             !  enddo
             !enddo
           else if(mraerosol) then
             do k=1,levs
               do i=1,im
                 dqdt_water_vapor(i,k)             = RQVBLTEN(i,k) !/(1.0 + qv(i,k))
                 dqdt_liquid_cloud(i,k)            = RQCBLTEN(i,k) !/(1.0 + qv(i,k))
                 dqdt_cloud_droplet_num_conc(i,k)  = RQNCBLTEN(i,k)
                 dqdt_ice_cloud(i,k)               = RQIBLTEN(i,k) !/(1.0 + qv(i,k))
                 dqdt_ice_num_conc(i,k)            = RQNIBLTEN(i,k)
               enddo
             enddo
             if(ldiag3d .and. .not. flag_for_pbl_generic_tend) then
               call dtend_helper(100+ntqv,RQVBLTEN)
               call dtend_helper(100+ntcw,RQCBLTEN)
               call dtend_helper(100+ntlnc,RQNCBLTEN)
               call dtend_helper(100+ntiw,RQIBLTEN)
               call dtend_helper(100+ntinc,RQNIBLTEN)
             endif
           else
             !Thompson (2008)
             do k=1,levs
               do i=1,im
                 dqdt_water_vapor(i,k)   = RQVBLTEN(i,k) !/(1.0 + qv(i,k))
                 dqdt_liquid_cloud(i,k)  = RQCBLTEN(i,k) !/(1.0 + qv(i,k))
                 dqdt_ice_cloud(i,k)     = RQIBLTEN(i,k) !/(1.0 + qv(i,k))
                 dqdt_ice_num_conc(i,k)  = RQNIBLTEN(i,k)
                 !dqdt_ozone(i,k)         = 0.0
               enddo
             enddo
             if(ldiag3d .and. .not. flag_for_pbl_generic_tend) then
               call dtend_helper(100+ntqv,RQVBLTEN)
               call dtend_helper(100+ntcw,RQCBLTEN)
               call dtend_helper(100+ntiw,RQIBLTEN)
               call dtend_helper(100+ntinc,RQNIBLTEN)
             endif
             !do k=1,levs
             !  do i=1,im
             !    qgrs_water_vapor(i,k)            = qgrs_water_vapor(i,k)    + (RQVBLTEN(i,k)/(1.0+RQVBLTEN(i,k)))*delt
             !    qgrs_liquid_cloud(i,k)           = qgrs_liquid_cloud(i,k)   + RQCBLTEN(i,k)*delt
             !    qgrs_ice_cloud(i,k)              = qgrs_ice_cloud(i,k)      + RQIBLTEN(i,k)*delt
             !    qgrs_cloud_ice_num_conc(i,k)     = qgrs_cloud_ice_num_conc(i,k)     + RQNIBLTEN(i,k)*delt
             !    !dqdt_ozone(i,k)        = 0.0
             !  enddo
             !enddo
           endif !end thompson choice
        elseif (imp_physics == imp_physics_nssl) then
           ! NSSL
             do k=1,levs
               do i=1,im
                 dqdt_water_vapor(i,k)             = RQVBLTEN(i,k) !/(1.0 + qv(i,k))
                 dqdt_liquid_cloud(i,k)            = RQCBLTEN(i,k) !/(1.0 + qv(i,k))
                 dqdt_cloud_droplet_num_conc(i,k)  = RQNCBLTEN(i,k)
                 dqdt_ice_cloud(i,k)               = RQIBLTEN(i,k) !/(1.0 + qv(i,k))
                 dqdt_ice_num_conc(i,k)            = RQNIBLTEN(i,k)
                 IF ( nssl_ccn_on ) THEN ! 
                   dqdt_cccn(i,k)      = RQNWFABLTEN(i,k)
                 ENDIF
               enddo
             enddo

        elseif (imp_physics == imp_physics_gfdl) then
           ! GFDL MP
           do k=1,levs
             do i=1,im
               dqdt_water_vapor(i,k)   = RQVBLTEN(i,k) !/(1.0 + qv(i,k))
               dqdt_liquid_cloud(i,k)  = RQCBLTEN(i,k) !/(1.0 + qv(i,k))
               dqdt_ice_cloud(i,k)     = RQIBLTEN(i,k) !/(1.0 + qv(i,k))
               !dqdt_rain(i,k)          = 0.0
               !dqdt_snow(i,k)          = 0.0
               !dqdt_graupel(i,k)       = 0.0
               !dqdt_ozone(i,k)         = 0.0
             enddo
           enddo
           if(ldiag3d .and. .not. flag_for_pbl_generic_tend) then
             call dtend_helper(100+ntqv,RQVBLTEN)
             call dtend_helper(100+ntcw,RQCBLTEN)
             call dtend_helper(100+ntiw,RQIBLTEN)
           endif
           !do k=1,levs
           !  do i=1,im
           !    qgrs_water_vapor(i,k)            = qgrs_water_vapor(i,k)    + (RQVBLTEN(i,k)/(1.0+RQVBLTEN(i,k)))*delt
           !    qgrs_liquid_cloud(i,k)           = qgrs_liquid_cloud(i,k)   + RQCBLTEN(i,k)*delt
           !    qgrs_ice_cloud(i,k)              = qgrs_ice_cloud(i,k)      + RQIBLTEN(i,k)*delt
           !    !dqdt_ozone(i,k)        = 0.0
           !  enddo
           !enddo
       else
!          print*,"In MYNN wrapper. Unknown microphysics scheme, imp_physics=",imp_physics
           do k=1,levs
             do i=1,im
               dqdt_water_vapor(i,k)   = RQVBLTEN(i,k) !/(1.0 + qv(i,k))
               dqdt_liquid_cloud(i,k)  = RQCBLTEN(i,k) !/(1.0 + qv(i,k))
               dqdt_ice_cloud(i,k)     = 0.0
               !dqdt_rain(i,k)          = 0.0
               !dqdt_snow(i,k)          = 0.0
               !dqdt_graupel(i,k)       = 0.0
               !dqdt_ozone(i,k)         = 0.0
             enddo
           enddo
           if(ldiag3d .and. .not. flag_for_pbl_generic_tend) then
             call dtend_helper(100+ntqv,RQVBLTEN)
             call dtend_helper(100+ntcw,RQCBLTEN)
             call dtend_helper(100+ntiw,RQIBLTEN)
           endif
       endif
       
       if (lprnt) then
          print*
          print*,"===Finished with mynn_bl_driver; output:"
          print*,"T:",t3d(1,1),t3d(1,2),t3d(1,levs)
          print*,"TH:",th(1,1),th(1,2),th(1,levs)
          print*,"rho:",rho(1,1),rho(1,2),rho(1,levs)
          print*,"exner:",exner(1,1),exner(1,2),exner(1,levs)
          print*,"prsl:",prsl(1,1),prsl(1,2),prsl(1,levs)
          print*,"dz:",dz(1,1),dz(1,2),dz(1,levs)
          print*,"u:",u(1,1),u(1,2),u(1,levs)
          print*,"v:",v(1,1),v(1,2),v(1,levs)
          print*,"sqv:",sqv(1,1),sqv(1,2),sqv(1,levs)
          print*,"sqc:",sqc(1,1),sqc(1,2),sqc(1,levs)
          print*,"sqi:",sqi(1,1),sqi(1,2),sqi(1,levs)
          print*,"rmol:",rmol(1)," ust:",ust(1)
          print*,"dx(1)=",dx(1),"initflag=",initflag
          print*,"Tsurf:",tsurf(1)," Thetasurf:",ts(1)
          print*,"HFX:",hfx(1)," qfx",qfx(1)
          print*,"qsfc:",qsfc(1)," ps:",ps(1)
          print*,"wspd:",wspd(1)," rb=",rb(1)
          print*,"znt:",znt(1)," delt=",delt
          print*,"im=",im," levs=",levs
          print*,"PBLH=",pblh(1)," KPBL=",KPBL(1)," xland=",xland(1)
          print*,"vdfg=",vdfg(1)," ch=",ch(1)
          !print*,"TKE:",TKE_PBL(1,1),TKE_PBL(1,2),TKE_PBL(1,levs)
          print*,"qke:",qke(1,1),qke(1,2),qke(1,levs)
          print*,"el_pbl:",el_pbl(1,1),el_pbl(1,2),el_pbl(1,levs)
          print*,"Sh3d:",Sh3d(1,1),sh3d(1,2),sh3d(1,levs)
          print*,"exch_h:",exch_h(1,1),exch_h(1,2),exch_h(1,levs)
          print*,"exch_m:",exch_m(1,1),exch_m(1,2),exch_m(1,levs)
          print*,"max cf_bl:",maxval(cldfra_bl(1,:))
          print*,"max qc_bl:",maxval(qc_bl(1,:))
          print*,"dtdt:",dtdt(1,1),dtdt(1,2),dtdt(1,levs)
          print*,"dudt:",dudt(1,1),dudt(1,2),dudt(1,levs)
          print*,"dvdt:",dvdt(1,1),dvdt(1,2),dvdt(1,levs)
          print*,"dqdt:",dqdt_water_vapor(1,1),dqdt_water_vapor(1,2),dqdt_water_vapor(1,levs)
          print*,"ktop_plume:",ktop_plume(1)," maxmf:",maxmf(1)
          print*,"nup:",nupdraft(1)
          print*
       endif

       if(allocated(save_qke_adv)) then
         if(ldiag3d .and. .not. flag_for_pbl_generic_tend) then
           idtend = dtidx(100+ntke,index_of_process_pbl)
           if(idtend>=1) then
             dtend(:,:,idtend) = dtend(:,:,idtend) + qke_adv-save_qke_adv
           endif
         endif
         deallocate(save_qke_adv)
       endif

  CONTAINS

    SUBROUTINE dtend_helper(itracer,field,mult)
      real(kind=kind_phys), intent(in) :: field(im,levs)
      real(kind=kind_phys), intent(in), optional :: mult(im,levs)
      integer, intent(in) :: itracer
      integer :: idtend
      
      idtend=dtidx(itracer,index_of_process_pbl)
      if(idtend>=1) then
        if(present(mult)) then
          dtend(:,:,idtend) = dtend(:,:,idtend) + field*dtf*mult
        else
          dtend(:,:,idtend) = dtend(:,:,idtend) + field*dtf
        endif
      endif
    END SUBROUTINE dtend_helper

! ==================================================================
  SUBROUTINE moisture_check2(kte, delt, dp, exner, &
                             qv, qc, qi, th        )
  !
  ! If qc < qcmin, qi < qimin, or qv < qvmin happens in any layer,
  ! force them to be larger than minimum value by (1) condensating 
  ! water vapor into liquid or ice, and (2) by transporting water vapor 
  ! from the very lower layer.
  ! 
  ! We then update the final state variables and tendencies associated
  ! with this correction. If any condensation happens, update theta/temperature too.
  ! Note that (qv,qc,qi,th) are the final state variables after
  ! applying corresponding input tendencies and corrective tendencies.

    implicit none
    integer,  intent(in)     :: kte
    real, intent(in)         :: delt
    real, dimension(kte), intent(in)     :: dp, exner
    real, dimension(kte), intent(inout)  :: qv, qc, qi, th
    integer   k
    real ::  dqc2, dqi2, dqv2, sum, aa, dum
    real, parameter :: qvmin1= 1e-8,    & !min at k=1
                       qvmin = 1e-20,   & !min above k=1
                       qcmin = 0.0,     &
                       qimin = 0.0

    do k = kte, 1, -1  ! From the top to the surface
       dqc2 = max(0.0, qcmin-qc(k)) !qc deficit (>=0)
       dqi2 = max(0.0, qimin-qi(k)) !qi deficit (>=0)

       !update species
       qc(k)  = qc(k)  +  dqc2
       qi(k)  = qi(k)  +  dqi2
       qv(k)  = qv(k)  -  dqc2 - dqi2
       !for theta
       !th(k)  = th(k)  +  xlvcp/exner(k)*dqc2 + &
       !                   xlscp/exner(k)*dqi2
       !for temperature
       th(k)  = th(k)  +  xlvcp*dqc2 + &
                          xlscp*dqi2

       !then fix qv if lending qv made it negative
       if (k .eq. 1) then
          dqv2   = max(0.0, qvmin1-qv(k)) !qv deficit (>=0)
          qv(k)  = qv(k)  + dqv2
          qv(k)  = max(qv(k),qvmin1)
          dqv2   = 0.0
       else
          dqv2   = max(0.0, qvmin-qv(k)) !qv deficit (>=0)
          qv(k)  = qv(k)  + dqv2
          qv(k-1)= qv(k-1)  - dqv2*dp(k)/dp(k-1)
          qv(k)  = max(qv(k),qvmin)
       endif
       qc(k) = max(qc(k),qcmin)
       qi(k) = max(qi(k),qimin)
    end do

    ! Extra moisture used to satisfy 'qv(1)>=qvmin' is proportionally
    ! extracted from all the layers that has 'qv > 2*qvmin'. This fully
    ! preserves column moisture.
    if( dqv2 .gt. 1.e-20 ) then
        sum = 0.0
        do k = 1, kte
           if( qv(k) .gt. 2.0*qvmin ) sum = sum + qv(k)*dp(k)
        enddo
        aa = dqv2*dp(1)/max(1.e-20,sum)
        if( aa .lt. 0.5 ) then
            do k = 1, kte
               if( qv(k) .gt. 2.0*qvmin ) then
                   dum    = aa*qv(k)
                   qv(k)  = qv(k) - dum
               endif
            enddo
        else
        ! For testing purposes only (not yet found in any output):
        !    write(*,*) 'Full moisture conservation is impossible'
        endif
    endif

    return

  END SUBROUTINE moisture_check2

  END SUBROUTINE mynnedmf_wrapper_run

!###=================================================================

END MODULE mynnedmf_wrapper
